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We show that the substitutional vacancy in graphene forms a dynamical Jahn- Teller center. 
The adiabatic potential surface resulting from the electron-lattice coupling was computed using 
density-functional methods and subsequently the Schrodinger equation was solved for the nuclear 
motion. Our calculations show a large tunneling splitting SF of about 86 cm^^. The effect results 
in a large delocalization of the carbon nuclear wave functions around the vacancy leading to a 
significant broadening of the Jahn- Teller active sp^a electron states. The tunneling splitting should 
be observable in electron paramagnetic resonance and two-photon resonance scattering experiments. 
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In spite of its deceptively simple honeycomb lattice 
structure, graphene has quickly become a new paradigm 
for testing a variety of ideas in condensed matter physics. 
The much celebrated linear band structure of graphen^ 
leads to a host of unusual behaviors such as Klein tun- 
neling, chiral electrons, minimum conductivity, negative 
refraction, half- integer quantum Hall effect, and new fea- 
tures in the Kondo and RKKY interactions leading to 
quantum criticality.l^Hs! Vacancies in the carbon based sys- 
tems have been of considerable interest for quite some 
time now, especially i n the c ontext of magnetism with- 
out magnetic atoms.l^ ^ l ^° l ^^ l Quite remarkably, it has 
been shown that a vacancy introduces a quasi-localized 
midgap state in the tt bands with ~ 1/r decay on account 
of the particle-hole symmetry.E^tSl An interesting conse- 
quence of this is the partial occupation of the vacancy- 
induced cr-band states, which leads then to a Jahn- Teller 
(JT) distortion. The JT distortion could be static or 
dynamic. In the latter, the potential barrier between 
the different equivalent minima in the nuclear configura- 
tion space is small enough that the nuclei tunnel between 
the various minima leading to several interesting effects, 
while in the static JT effect, the nuclei are stuck to one 
minima or the other. In this Letter, we show that the 
vacancy forms a dynamical JT center in graphene ow- 
ing to the small quantum mechanical barrier for nuclear 
tunneling. 

Density-functional calculations^'^' show that the va- 
cancy introduces four electrons into the graphene bands 
as illustrated in Fig. [l] The JT effect comes from the 
partial occupation of the doubly-degenerate sp^a dan- 
gling bond states on the carbon triangle surrounding the 
vacancy and their coupling to the two vibrational modes 
of the triangle, given by the E ® e JT Hamiltoniar l^^ * ^^ ! 

+ g (QA + + G [{Ql - Qj) f, + 2Q1Q2T,] , 

where the terms are the nuclear kinetic energy, the elas- 
tic energy, and the linear and the quadratic JT cou- 
pling terms. Here the pseudospin r describes the two 
JT active, doubly-degenerate electronic states and 




FIG. 1: (Color online) Vacancy induced a and tt electron 
states (Vct and V,r) with the occupied states shown by dots 
with arrows (left). The nominal 2fiB {S = 1) magnetic mo- 
ment due to the localized states is reduced substantially due 
to the anti-ferromagnetic spin polarization of the band states, 
indicated by tt; tii in the local neighborhood of the vacancy. 
Right part shows the JT active electron states, \vi) and \v2), 
and the vibrational modes of the carbon triangle that they 
couple to. (Ti denotes the dangling sp^a bond orbital on a 
carbon atom adjacent to the vacancy. 



\v2) originating from the three sp^a dangling bonds 
on the carbon triangle: \vo) = (cri + CT2 + cr3)/"\/3, 

\vi) = + ^3)7^2, \V2) = (2cti - CT2 - (J3)/V6, 

with energies Eq — —2t and £'1^2 = t and symme- 
tries Ai and E, respectively, with the —t being the 
cr-electron hopping between the neighboring sites on 
the triangle, and \vi) transforms like x and \v2) like 
y. On the other hand, the Pz orbitals, responsible for 
the linear 'tt' Dirac bands, introduce the quasi-localized 
midgap state, which becomes singly-occupied due to 
Hund's coupling, leaving a lone electron to occupy the 
CT-derived doubly-degenerate E state. This explains the 
the relative positions of the vacancy states shown in 
Fig. [T] Turning now to the three vibrational modes 
of the triangle: |Qo) = (0, 2, v^, -1, -^3, -l)/yT2, 
IQi) = (0,2,-73,-1,73,-1)/^^, IQ2) = (2,0,-1, 
a/3, —1, -'V^)/\/l2^Qo is the stretching mode and the 
doubly-degenerate Qi and Q2 modes are JT active, split- 



2 



ting the upper two Vcr bands as shown in Fig. [T] The pa- 
rameters in the Hamiltonian are the carbon mass A/, the 
elastic energy and the hnear and quadratic JT cou- 
pHng parameters g and G, respectively. Diagonalization 
of the potential terms in Eq. [T] leads to the well-known 
adiabatic potential surface (APS) for the nuclear motion 

E± = ^Kp^ ± py/g^ + G^p^ + 2gGpcos{3<P), (2) 

where p — \JQ\+Q\ and — tan^^(Q2/Qi) are the po- 
lar coordinates and E± denote the two potential sheets. 
Without the quadratic coupling (G = 0), one gets the 
Mexican hat APS, while with it we have three minima in 
the (Qi, Q2) space (Fig- [2]). The electronic eigenfunction 
for the lower sheet ia^^ 

W) ^ [sin(0/2)|i;i) + cos((/./2)|t-2)] x e*^/^, (3) 

where the phase factor assures single-valuedness as one 
moves around the origin and leads to a Berry phase. 

In order to study the APS, we have computed the to- 
tal energy as a function of the vibronic coordinates using 
the spin-polarized density functional all-electron linear 
augmented plane- waves (LAPW) methodi^ and the gra- 
dient approximation (GGA) for the exchange-correlation 
functional .fi^ We used a 32-atom supercell with a single 
vacancy and obtained a fully relaxed structure, which 
yielded a planar structure with an isosceles triangle for 
the carbon atoms surrounding the vacancy with two long 
bonds (2.66 A) and one short bond (2.41 A). This is 
equivalent to the distortion: Qo = 0.08 A, Qi = 0.166 
A, and Q2 — 0. We then took a series of structures 
with varying distortions, Qi and Q2, and in each case 
optimized the rest of the carbon atoms in the supercell. 
We note that while the literature is divided regarding 
whether the relaxed structure with a vacancy is planar 
or non-planar, the three- fold symmetry of the adiabatic 
potential surface occurs in either case, being tied to the 
symmetry of the honeycomb lattice itself. The calculated 
energies are shown in Fig. [2] which yields the JT distor- 
tion radius po — 0.165 A, the JT stabilization energy 
iJjT= 110 meV and the tunneling barrier height /3=19 
meV. Comparison of these results with Eq. ([2| yields 
the stiffness constant K = 9.3 eV j h? and the linear 
and the quadratic JT parameters g = 1.46 eV j A and 
G = 0.38 eV I A^, respectively. For the case of LaMnOs, 
a well-known system with a strong JT interaction, while 
the K and q are about the same, the warping parameter 
G = 2.0 eVj A2 is significantly large ,'22' which results in 
a static JT effect with the nuclei stuck to one potential 
minimum. In contrast, the weaker warping term G in 
graphene leads to a small barrier height for nuclear tun- 
neling and consequently to the dynamic JT effect, where 
the nuclei tunnel between the three minima in the APS. 
Since the phonon frequency for the nuclear motion in the 
potential well hu) « 57 meV is much larger than the 
barrier /3, the nuclei cannot be localized in one of the 
potential wells. 
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FIG. 2: (Color online) Total energy as a function of the vi- 
bronic distortion Q\ computed from the DFT (red dots) and 
fitted to the adiabatic energy E- in Eq. (|2|(full line) (Top). 
The triangles indicate the configurations at the three extrema. 
Bottom figure shows the corresponding energy contours in the 
Qi — Q2 plane, with the three equivalent minima (dots) sepa- 
rated by the tunneling barriers (crosses). The contour values 
are: -0.11 -I- 0.001 x (2") in units of eV, where n = 0, 1, 7 
labels the contours and F denotes the nuclear hopping integral 
in the tight binding description. 



It is difficult to treat the dynamical JT effect using 
DFT when many vibrational modes are present as in case 
of a JT center in the crystal and often the sing le-mo de 
approximation is made with remarkable success. 'i^'^^ In 
the present case, due to the localized nature of the JT- 
active states (dangling sp^ bond orbitals pointed towards 
the vacancy) , the JT coupling to modes belonging to fur- 
ther neighbor shells is weak (for the second shell, we find 
g' ~ g/&) and since the higher shell stiffness constants are 
significantly larger than for the first shell for the vacancy 
center, the single-mode approximation captures the es- 
sential physics in the present case. 

The basic features of the collective nuclear-electronic 
motion may be described by adopting a simple tight- 
binding approach, familiar from the electronic struc- 
ture theory. We write the collective wave function 
as the linear combination I^P) = ci (j)i{R) ipl{R,r) + 
C2 MR) V'|(^,?-) + C3 MR) i'liR^r), where R{r) is the 
nuclear (electronic) coordinate, (f>i{R) solves the nuclear 
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Schrodinger equation in the vicinity of the potential min- 
ima, 



[Tr + V,{R)]UR) = Eo(t>,{R), 



(4) 



and ipi{R, r) satisfies the electronic Schrodinger equation 
for the fixed nuclear position R = (Qi,Q2)- The elec- 
tronic wave function is restricted to the Hilbert space 
{\vi),\v2)) and has the form Eq. ^ for a given nu- 
clear coordinate R. Thus the energy eigenstates assume 
the Born-Oppenheimer form \'^{R)) = <l'„(i?)|'0e(i?, r)), 
where = ci0i(_R) -t- C2<j)2{R) + cz4'3,{R) is a hnear 

combination of the nuclear orbitals. The eigenstates can 
then be obtained from the diagonalization of the 3x3 
Hamiltonian 



H 




(5) 



where the phase factor e*"^ will be discussed momen- 
tarily, Eq is the on-site energy, and T is the nu- 
clear hopping integral in the adiabatic approximation 
r = {ct)^{R)r{R.r)\/W{KM2{R)r{R,r)) « -AV^ x 
F. In obtaining the last result, the normaliza- 
tion {ip'^ {R, r)\ip'^ {R, r)) = 1 has been used, F — 
J (f>l{R)(f>2{R)cfiR is the Frank-Condon factor, and the 
deviation of the lower APS potential from the well po- 
tential, AV{R) = V-^{R) — Vi{R), has been approximated 
by its value —AV at the saddle point (marked by a cross 
in the bottom panel in Fig ([2])), since that's where most 
of the contribution to the integral comes from. 

The magnitude of the nuclear hopping F may be esti- 
mated by assuming a one-dimensional motion of the nu- 
clei in the azimuthal direction, along the circle of radius 
Pq and by computing the quantities AV and F. The ID 
motion is reasonable since by expanding the adiabatic po- 
tential V- around the potential minima, the spring con- 
stant for azimuthal motion is found to be K' = 9G, which 
is about half of the spring constant K for radial motion. 
This corresponds to a phonon frequency of hiu « 58 meV 
for radial motion and « 34 meV for the azimuthal mo- 
tion. The latter is of the same order of magnitude as 
the tunneling barrier of 19 meV, which again indicates 
strong tunneling between the three minima. Now, taking 
the nuclear wave functions as the ID simple harmonic os- 
cillator wave function localized at the potential minima: 
cl){x) = (7rZ2)-i/4exp [-x^/i2P)], where / = H/VMK' 
and X is the length along the azimuthal direction, the 
Frank-Condon factor becomes simply the overlap integral 
between two displaced harmonic oscillator wave func- 
tions, with the result: F = 2~^/^exp [— a^/(4/^)], where 
a — I'KpQji is the distance between two minima along 
the circle. Meanwhile, the potential difference between 
the minimum and the saddle point can be found to be 
AV — [tt'^K' /1% — 2G). Plugging in the numerical 
values, we find F ~ 0.13 and AV^ ~ 0.035 eV, so that the 
hopping integral F « AV^ x F ~ —37 cm~^. 

Finally, in addition to the hopping integral, the adia- 
batic motion of the electron results in a fictitious mag- 



netic field seen by the nuclei with the vector potentiaPi' 

h 



A 



q 



lm{iP,{R, r)|Vi?^e(i?,r)), 



(6) 



which adds a phase factor to the hopping amplitude in 
the Hamiltonian ([5|. The modified hopping in the pres- 
ence of the magnetic field, from point a to b, is given by 
the expressioip2l 



(7) 



It immediately follows from Eqs. (|3| and (|6|) that A = 
—2^^hq~^e^, so that the phase factor in Eq. (l7| is sim- 



ply e 



pi7r/3 



This phase factor is very important as 



without this, the symmetry of the ground state is in- 
correctly predicted. Diagonalization of the Hamiltonian 
Eq. ^ with the correct phase factor yields a doubly- 
degenerate nuclear ground-state with energy F, with the 
singly-degenerate excited state at energy — 2F, so that 
the energy separation is 3|F| = 111 cm^^. 
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FIG. 3: (Color online) Eigenvalues obtained by diagonaliza- 
tion of Eq. ([T]) using the basis set Eq. ([8| as a function of 
the scaled coupling strengths Xg and AG. Numbers inside the 
figure indicates the degeneracies. For A = 0, eigenstates of 
the two-dimensional harmonic oscillator are reproduced. 

This crude but conceptually rich tight-binding result 
may be compared to the exact, brute-force diagonaliza- 
tion of the full Hamiltonian Eq. ([!]) by expanding the 
combined nuc lear-e lectronic wave function j^I') in a com- 
plete basis selP^'^Il 
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A 

i.c\r (4) 



(4)" iAr 



\vi) 



(8) 



where c\ , create harmonic oscillator states along the 
Q11Q2 axes centered at the origin and Anm and Bnm 
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are the expansion coefficients. This procedure requires 
no additional consideration of a fictitious magnetic field 
and also yields the full solutions in addition to the lowest 
three states obtained from the tight-binding theory. The 
results are shown in Fig. ([s]). The magnitude of the 
tunneling splitting 3|r| — 86 cm^^ compares very well 
with the tight-binding result. 

The large value of the tunneling splitting as compared 
to the strain splitting, typical value of whiclP^ is (5 ~ 10 
cm~^, results in the delocalization of the nuclear wave 
function. If the reverse were true, then the nuclei would 
be more or less stuck in one or the other potential well due 
to the removal of the degeneracy of the three APS minima 
by the local strain caused by the invariable presence of 
defects. This would therefore lead to a static distortion 
of the nuclear framework resulting in the static JT effect. 
For the dynamical JT effect, the tunneling splitting must 
be strong enough to overcome the strain splitting, so that 
the nuclei can tunnel between all APS minima, which is 
the case for graphene. 

The nuclear probability density in real space |^'jv('')P 
can be computed from the corresponding quantity in the 
configuration space 

(9) 

where 0„ is the vl^ harmonic oscillator eigenfunction. 
The calculated |5'Ar(r)p is shown in Fig. Q, which 
indicates a significant spread of the nuclear wave func- 
tion of the carbon triangle, about 0.1 A from the equi- 
librium positions. We note that this is not washed out 
by the lattice thermal vibrations, which causes the nu- 
clear vibrational amplitude, estimated from the expres- 
sion \KQ^ — ^ksT to be about 0.05 A at room temper- 
ature. 

The spread of the nuclear wave function broadens the 
energy of the JT split electronic states as well, so that 
they are not sharp (5-function states any longer. In 
the adiabatic approximation, the electronic density-of- 
states is given by p{E) = J^q.q^ I*w(<3i, <32)P x [S{E- 
S-{Qi,Q2)) + S{E — e+{Qi,Q2)), where e± denote the 
two JT split states as in the expression ^ without the 
elastic energy term. The results are shown in the inset 
of Fig. Q, which predicts a rather large width, of the 
order of 0.15 eV, due to the JT effect. Thus these states 
should appear as rather broad states in scanning tunnel- 
ing experiments. In contrast to this, the broadening of 
the midgap Vtt state is expected to be rather small. In 
fact, it is exactly zero if only the nearest-neighbor hop- 
ping is retained. This is borne out by the less than 
5 meV width of the midgap state, seen in the scanning 
tunneling experiments.'^^ 

In conclusion, we showed that the substitutional va- 
cancy in graphene forms a dynamical JT center owing to 
a weak potential barrier for tunneling between the three 
minima in the adiabatic potential surface. The doubly- 
degenerate nuclear ground state with the tunneling split- 
ting of about 86 cm~^ originates from the combined 




FIG. 4: (Color online) Nuclear probability density |^'jv(r)p 
showing the symmetric distortion of the carbon atoms from 
the ideal position of an equilateral triangle (solid line). The 
nuclei move in a correlated manner so that the most probable 
configuration is one of the three isosceles triangles (dashed 
lines) corresponding to the three minima of the APS. The nu- 
clear motion of the nearby atoms show much smaller devia- 
tion from their equilibrium positions. Inset shows a significant 
broadening, computed within the adiabatic approximation, of 
the JT active electron states due to the spread of the nuclear 
wave function. 



nuclear-electronic motion, which may be cast in terms of 
a Berry phase acquired due to a fictitious magnetic field 
experienced by the nuclei caused by the adiabatic mo- 
tion of the electrons. The splitting should be observable 
in the electron paramagnetic resonance and two-photon 
resonance scattering experiments, which have been used 
to study the JT effects in the triatomic molecules. The 
quantum mechanical spread of the nuclear wave func- 
tion is predicted to lead to a significant broadening of 
the JT split dangling bond states. Recently, it has been 
proposecP^ that the entanglement between the nuclear 
and electronic motion in a dynamical JT system may be 
exploited in quantum computation, leading to the possi- 
bility of yet another novel application for graphene. 
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